Chapter 4 

Methods of Statistical Mechanics 



4.1 Mean-Field Theories 

Very few models of statistical mechanics have been solved exactly; in most of 
the cases one has to rely on approximative methods. Among them, the mean-field 
approximation (MFA) is one of the most widely used. The advantage of the mean- 
field theory is its simplicity and that it correctly predicts the qualitative features of 
a system in most cases. 

The essence of the mean-field theory is the assumption of statistical indepen- 
dence of the local ordering (spins in the case of magnetic systems). The interac- 
tion terms in the Hamiltonian are replaced by an effective, "mean field" term. In 
this way, all the information on correlations in the fluctuations is lost. Therefore, 
the mean-field theory is usually inadequate in the critical regime. It usually gives 
wrong critical exponents, in particular at low spatial dimensions when the number 
of NN is small. 

The MFA becomes exact in the limit as the number of interacting neighbours 
z — > oo. This is the case when the range of interaction z — > oo or if the spatial 
dimension is high. In both cases, a site "feels" contributions from many neigh- 
bours and therefore the fluctuations average out or become even irrelevant (for 
Ising systems, they become irrelevant in d > 4). At low d, the MFA must be used 
with caution. Not only that it predicts wrong critical exponents, its predictions are 
even qualitatively wrong sometimes. For example, it predicts long-range order 
and finite critical temperature for the d = 1 Ising model. In the following, we 
shall introduce and formulate the MFA for magnetic systems. In the Appendix A 
we summarize basic thermodynamic relations for magnetic systems. The results, 
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however, are valid also for other systems. 

Two formulations of the mean-field approximation. We will now formulate 
the mean-field approximation in two different ways. 

• A: The "Weiss molecular field theory" according to its "inventor" Pierre 
Weiss. This method is straightforward and easy to understand. It gives us 
an expression for the order parameter but does not tell us anything about the 
free energy or partition function. 

• B: The second method will be the Bragg-Williams approximation which is 
more sophisticated and is based on the free energy minimization. 

We will develop the mean-field approximation on the Ising model in an exter- 
nal magnetic field. 

4.1.1 Weiss Molecular Field of an Ising System 

Weiss formulated a theory of ferromagnetism in which he assumed that the effect 
of the neighbouring spins on a given spin can be described by a fictitious "molec- 
ular field" which is proportional to the average magnetization of the system. 
We start with the Ising Hamiltonian in an external magnetic field: 

(ij) i 

The local field acting on the spin at the site % is: 

Hi = Jj2Sj + H - (4-2) 

j 

It depends on the configuration of all the neighbouring spins around the site i. 
If the lattice is such that each spin has many neighbours, then it is not a bad 
approximation to replace the actual value of the neighbouring spins by the mean 
value of all spins, 

Y^S^zm, (4.3) 

3 

where z is the number of nearest neighbours and m is the mean (average) magne- 
tization per site, m = ~^J2j Sj- With (4.3), the local field becomes 

Hi = zJm + H, (4.4) 
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Figure 4.1: Temperature dependence of the order parameter in the mean-field 
approximation. 

The actual local field acting on the site % was replaced by its mean value and is 
independent of the site, Hi — > H M f- In the mean-field approximation the original 
Hamiltonian is replaced by a mean-field Hamiltonian 



which is equivalent to the Hamiltonian of a non-interacting spin system in an 
external magnetic field and can be solved exactly: 



This is a self-consistent, transcendental equation for m. It has a non-trivial solu- 
tion (m 7^ 0) for (3z J > 1, i.e., at low temperatures. The temperature dependence 
of the order parameter m for H = is shown in Fig. 4. 1. It vanishes at the critical 
temperature Tc = zJ. 

4.1.2 Bragg- Williams Approximation 

The Gibbs free energy is a function of T and H, 



7~Lmf — — 




(4.5) 



m = tanh(/5 H M f) = tanh[/3(z Jm + H)]. 



(4.6) 



dG = -SdT - MdH. 



(4.7) 
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This equation tells us that for fixed T and H, 5G = 0. In equilibrium, thus, G 
is an extremum of a trial free energy Q with respect to any variational parameter. 
We will now search the trial free energy Q. 

We start with the density matrix and calculate the Gibbs free energy of an 
Ising system. In the Bragg- Williams (mean-field) approach we approximate the 
total density matrix p by a direct product of the density matrices for individual 
spins, 

P^UpitxUe-^. (4.8) 

i i 

Physically this means that the spin fluctuations are considered to be uncorrelated, 
the orientation of a spin is statistically independent of the orientation of any other 
spin. There are no correlations in the spin orientation except those coming from 
the long-range order. Under this assumption, the Ising Hamiltonian is not only 
diagonal in the spin space (the state Si = +1 is uncoupled from the state Si = 
— 1; physically this means that the Ising system has no dynamics), but it is also 
independent of all other spins in the system. Thus, the trial density matrices p, 
are diagonal and of the general form (remember, its trace must be = 1): 



Ul + m') 



2 



\{l -m'i) 



(4.9) 



They depend on the variational parameter which is site-independent, m! i = ml, 
because of translational invariance. With this density matrix the average spin 
orientation is: 

{S i )=Tr(S iPi )=m. (4.10) 

Because of statistical independence of the spins in MFA, the short-range spin 
correlations are: 

(S i S j ) = (S i )(S j )=m 2 . (4.11) 
We can write the (trial) enthalpy E = (H) [where H is given by the Eq. (4.1)] as: 

E = -^zJNm' 2 - NHm' (4.12) 

and the trial entropy as: 

, m t i \ at, fl + m' 1 + m' l—m' 1 - m'\ 
S = -k B Tr(p In p) = -Nk B — - — In — 1 — In — - — 
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Figure 4.2: Trial Gibbs free energy g(T,H;m') vs. w! at H = 0.1 z J. Solid 
circles denote stable states, open circle a metastable state, and the solid triangle an 
unstable state. The curves are labelled by T/T c . m and G(T, H) are determined 
by the minimum of Q. 

(4.13) 

With these expressions, the trial Gibbs free energy becomes: 

g(H,T;m') = E-TS = E + k B T Tr(plnp) = -\zJNm' 2 - NHm'+ 

V m' l + m' 1 — m' 1 — m' \ A A 

In 1 In (4.14) 

2 2 2 2 ; v ' 

It depends on H, T, and on the variational parameter ml. Q for H = 0.1 (in 
units of z J) is shown in Fig. 4.2. To get the true equilibrium Gibbs free energy, 
G(H, T; m) has to be minimized with respect to the variational parameter ml, 

G(T,H) = min Q(H, T;m'). (4.15) 

m' 

The minimization gives: 

j 1 _|_ rri 

zJm + H = -fc^Tln — — — , (4.16) 
2 1 — m 
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Figure 4.3: External-field dependence of the Ising-model magnetization in the 
mean-field approximation. T c is the critical temperature. Solid line: stable states, 
dashed line: metastable states, dotted line: unstable states. 



This is the familiar "equation of state" of magnetic systems. 

The field-dependence of m, Eq. (4.17) is shown in Fig. 4.3. At high tem- 
peratures (T > T c ), the order parameter m behaves as in a paramagnet, this is 
the paramagnetic phase. The spins are partially ordered only because of H. For 
H = 0, the spins are disordered, therefore this phase is also called the disordered 
phase. At low temperature (T < T c ), the spins are spontaneously ordered even 
in the absence of H. This is the spontaneously ordered phase. The system has 
long-range order for J > 0, all the spins are oriented predominantly in the same 
direction, this is the ferromagnetic phase. 

The Gibbs free energy is: 



or: 



m = tanh [f3(zJm + H)\. 



(4.17) 



G(T, H) 
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where m is related to H via the Eq. (4.17). The equation of state is 

(dG(T, H)\ 



M = 



dH 



(M = Nm). 



(4.19) 



The Helmholtz free energy is F(T, M) = G(T, H) + EM: 



F MF (T, M) = Nk B T 



1 + m 1 + m 
In 



1 — m , 1 — m 
H In 



1 



NzJm 2 . 



(4.20) 



The field dependence of the equilibrium Gibbs free energy G is shown in Fig. 
4.4. The equation of state now reads 



H 



fdF(T, MY 
\ dM 



(4.21) 



The Helmholtz free energy F(T, M, N) is shown in Fig. 4.5. 

As before, the critical temperature, where the spontaneous ordering vanishes 
in the absence of external magnetic field H as T is increased, is determined by the 
condition f3 c zJ = 1, 



k B T c = z\J\ 



(4.22) 



In the past, many improvements to the simple MF theory presented here have 
been proposed. With these improvements, the MF critical temperature came closer 
to the exact critical temperature. None of these theories, however, was able to 
bring an improvements to the critical exponents, they kept their mean-field values. 
The real break-through brought the renormalization-group theory, which will be 
introduced in Section 4.3. 



4.1.3 Landau Theory of Continuous Phase Transitions 

In the previous Section we have seen that the trial Gibbs free energy Q(T, H; M') 
was an analytic function which had approximately a parabolic shape as a function 
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Figure 4.4: Equilibrium Gibbs free energy. At low temperatures, when a first- 
order transition at H = takes place, part of the curve corresponds to metastable 
states and part to unstable states. 

of the variational parameter M' above the critical temperature and a double min- 
imum below Tq. Close to the critical point and for H = 0, the minimum of Q is 
close to ml — and we can expand the trial mean-field free energy Q to a very 
good approximation in a power series of m. 



At T > T c , r must be > and below T c , it must be < 0, thus r changes sign and 
in most of the cases it can be approximated by a linear temperature dependence, 



The Landau theory is phenomenological and deals only with macroscopic 
quantities. The theory does not deal correctly with the fluctuations and is a mean- 
field theory. In fact, the Landau theory of continuous phase transitions is much 
more than just simple expansion of Q. It is based on the symmetry of the system. 
Landau postulated that the free energy can be expanded as a power series in m 
where only those terms that are compatible with the symmetry of the system are 
allowed. 




(4.23) 



r = r'(T-T c )/T c . 
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Figure 4.5: Equilibrium Helmholtz free energy. At low temperatures, where the 
transition is first order, F has two minima, corresponding to ±m. The region 
between the two minima is the coexistence region when the system can be partially 
in the +m, and partially in the — m state (domain structure!). 
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Formulation of the Landau Theory 

Critical anomalies take place on long lenthscales. Therefore it is enough and in 
order to change from discrete to continuous description: 

rrii — * m{f) 
hi ->• h(r). 

In formulating the Landau theory let us consider a macroscopic system and 
assume that there exists a (generalized; not necessarily magnetic) order parameter 
density m(r) averaged over microscopic, atomic distances a. Because of the aver- 
age, m(r) is now a continuous function of r (continuous in m and in r)\ There can 
be no fluctuations in m(f) on the lengthscale smaller than a (a^ lattice spacing). 
Correspondingly, m(r) has no Fourier components with wavevector greater than 
k = 2n/a. 

In lattice models, the canonical partition function appropriate to the Gibbs free 
energy of a spin system is: 

Y (T, H) = (4.24) 

{Si] 

(we have to sum over all spin configurations). Now, since m(r) is a continuous 
function of the coordinate r, we must replace the sum by a functional integral over 
all functions m(r). When doing this replacement one must not forget that m(r) 
is equal to the average of a microscopic spin Si over a small volume. That means 
that there can be several microscopic spin configurations which average into the 
same m(r). Let there be Q(m(r)) spin configurations that lead to the same m(r). 
Then, the sum over all spin configurations is equal to the functional integral over 
all functions m(r) where each function m(r) is weighted by Q(m(r)): 

]T -> f Dm(r) Q(m(r)). (4.25) 

The partition function is: 

Y(T,H) = J Dm(r)0(m(r)) e -^w(H;m(r))_ (4 26) 

Using the Boltzmann formula S = ks In Q, we can write Y in the form: 

Y(T,H) = jDmtfe-W&MryTS) 
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= /" Dm(r) e-^ T>JT;m < r » (4.27) 
where we have introduced the Landau free energy: 



Q(T, H; m(r)) = H(H; m(r) - TS(m(r)). 



(4.28) 



Notice the analogy between this equation and the trial Gibbs free energy ob- 
tained previously with the Bragg-Williams method! The Landau free energy is in 
fact a generalization of the trial Gibbs free energy. 

Notice also that the entropy S(m(r)) is not equal to the usual microcanonical 
entropy S (related to the number of configurations in an energy interval AE) 
because we did not make any summation over the energies! S(m(r)) is related 
to the number of microscopic configurations leading to the same function m{r). 
Only after integration over all the states in a certain energy interval we would 
get the usual thermodynamic entropy S. Therefore we can call S(m(r)) a partial 
entropy. 

We still have to work out the functional integral in order to get the (equilib- 
rium) partition function Y(T, H) ! From Y(T, H) the Gibbs free energy is: 

G(T, H) = -k B T In Y(T, H) = -k B T In J Dm e -^(T,ff;m)_ (4.29) 

In the Landau approach we write down the Landau free energy in a most gen- 
eral form as an expansion in terms of the order parameter and its derivatives by 
taking into account all the symmetry-allowed terms up to a given order: 

Q(T, H,m) — J dV (^-rm 2 + um 4 + ■ ■ ■ — Hm 

+ -#|Vm| 2 + ---) . (4.30) 

The highest-order coefficient in the expansion must always be positive, else the 
system were unstable. For the Ising model, this form can also be obtained by 
expanding G, Eq. (4.14), in powers of m. The last term appears if m is site- 
dependent: 



SiSj — > m(r) m(r + a) = m(r) m(r) + m(r) aVm(r) H — a 2 V 2 m(r) + 

2 



(4.31) 
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The first term is already included in the first line of (4.30), the second term van- 
ishes in systems with center of inversion (symmetry argument!) and the last term, 
after partial integration over the whole volume, gives the term proportional to 
| Vm| 2 . The prefactor g must also be positive always. A negative g would lead to 
completely unphysical wild fluctuations in m on extremely short lengthscales. 

Here we have assumed a one-component order parameter. In general, the order 
parameter can have several components or be complex. The coefficients are in 
principle temperature-dependent and depend also on the symmetry of the system. 
For example, the coefficient of the cubic term vanishes if the system is invariant 
under the reversal of m(r) when H = 0. Usually we assume that u is temperature 
independent whereas the quadratic term has a linear temperature dependence, 

r = r't, t = {T-T c )/T c (4.32) 

In the Landau theory we assume that the fluctuations of the order parameter 
are small, that the important values of m lie in a narrow range near its equilibrium 
value that minimizes the Landau free energy. In fact, the maximal contribution 
to the partition function comes from the term with minimal Landau free energy. 
In the thermodynamic limit, the probability of any other configuration will vanish 
rapidly. Therefore we make use of the saddle-point approximation, the functional 
integral (4.29) is equal to the maximum value of the integrand, multiplied by a 
constant. For constant field H, the functional form of m(r) that maximizes the 
integrand (minimizes the Landau free energy), is a constant, determined by the 
conditions 

dg d 2 g 

am om z 

In the saddle-point approximation we immediately obtain G(T, H) from Eq. (4.29): 

G(T, H) = G + min m g(T, H; m) = (4.34) 

In equilibrium, the system is at the bottom of the Landau free energy g. This 
determines m and G(T, H). From G(T, H) we get other thermodynamic quanti- 
ties, as before. 

The advantage of the Landau theory is that it allows one to study a macro- 
scopic system without knowing all its microscopic details. It is based on the sym- 
metry properties of the macroscopic system alone. 
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4.1.4 Critical Exponents in the Landau Theory 

We will now calculate the critical exponents of the Landau model. To be more 
specific, let us investigate the Landau model of magnetic phase transitions, when 

g = G + Qrm' 2 + urn' 4 - m'Hj V (4.35) 

In equilibrium, when dQ /dm! = 0, we get the equation of state: 

rm + 4um 3 - H = (4.36) 

• Order-parameter exponent f3 ( H — 0). For H = 0, and r > 0, the equation 
of state has a single solution, m — 0, there is no spontaneous magnetization 
above T c , as expected. For t < 0, the equation of state has three roots when 
H = 0. The root at m = belongs to an unstable state, so we are left with 
two real roots. The result is: 



r 



m = ±J^-t 1/2 j3 =\ t<0 



k m = £ > 

• Critical isotherm exponent 5 (T — T c ). Now, r = 0, and we have 4wm 3 = 
H, from which it follows 

m=\— =>- 5 = 3 £ = (4.38) 

V 4-u 

• Susceptibility exponent 7. The isothermal susceptibility is defined as 

(dm\ 

From the equation of state, Eq. (4.36), we get: 

r (wj T +12 ""* {m) t - 1 = (4 - 40) 

from which we find: 

— oc r 1 =^ 7 = 1 £ > 

f r'£ 

XT = I , (4.41) 

oc |£| _1 =>-7 = l £<0 
-2r'£ 1 1 
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Specific-heat exponent a (H — 0). The specific heat at constant (zero) field 
is 

<*-t(£) b (4.42, 

For t > 0, m = 0, and G = G = const, (that means that G does not 
depend on the reduced temperature t), so ch is zero. For t < (and H = 0), 

m 2 = — r/Au and 



„/2 

16m 



G = G - —t 2 (4.43) 



It follows: 



=>. a = t > 



c ff = r /2 T c „ „ (4.44) 



8m 



a = t < 



Comment 1: We see that the Landau theory gives mean-field (also called "classi- 
cal") critical exponents. 

Comment 2: At H ^ 0, there is no phase transition, all the quantities vary 
smoothly, only the exponent 5 is defined. 



4.1.5 Short-Range Order and Fluctuations 

Although the correlations and fluctuations are neglected in the MFA, we can still 
get an estimate of the short-range correlations. To obtain the mean-field values of 
the exponents v and r\, which describe the behaviour of the correlation function, 
we have to include the gradient term in the Landau free energy. In the frame- 
work of the mean-field theory, this leads to the Ornstein-Zernike approximation. 
To investigate the correlations, we consider the response of m(r) to a perturbing 
field acting, say, at the origin, H(r) = \8{r) via the "response function", the 
susceptibility: 

m(f) = J d d r'x(f-f)H(f). (4.45) 
The Fourier transform of the last convolution is simply: 

m(q) = x(0)H(q). (4.46) 
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To the lowest order in m and Vm, the Landau free energy is 
g(H;m(r)) = G + J d 



-g\Vm{r)\ 2 + -rm 2 (r) 
+um 4 (r) — Xm(r)S(r) 



(4.47) 



We first make Fourier transforms to the d— dimensional momentum space, 

m(q) = J d d r e _< * r m(f) 



m(r) 



(4.48) 



The Fourier transform of Vm is: 



Vm(f) = 



(2tt) c 



d d q (ye^)m{q) 
^— d J d d q iqe^ r m(q) 



(27T) 

The Fourier transform of a 5— function is a constant 

In the momentum space, the Landau free energy is: 
g(H;m(q))=G + J d d q 



(4.49) 



(4.50) 



1 1 

-gq 2 m 2 (q) + -rm 2 (q) 

+um 4 (q) - Xm(q) 



(4.51) 



Here, we approximated the Fourier transform of m 4 (f) by m 4 (g) ! 

In the mean-field approximation the modes with different momenta q are as- 
sumed to be independent of each other. Therefore we can minimize Q for each 
individual mode separately, 



dQ 



0, 



dm(q) 

and the equation of state in momentum space is: 

(gq 2 + r)m(q) + Aum 3 (q) = A. 



(4.52) 



(4.53) 
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For t > 0, when m vanishes in the absence of the field, we can neglect the cubic 
term, and we get: 

m(q) = —2^—, X(q) = (4-54) 

gqz _|_ r ggz _|_ r 

After a (somewhat tedious) Fourier transform back to real space we find in the 
limit of large r : 

x(r) oc m(r) oc r 2 ~ d exp(— r/£) =>- r\ = (4.55) 

where 

f = \/? ^ = i (4-56) 

Below 7b, m ^ 0, we replace m(f) by m+5m(r), where m is the spontaneous 
magnetization and 8m the response to the infinitesimal perturbing field A. We 
expand in powers of 5m and retain only linear-order terms. We find: 

Sm(q) = / - (4.57) 

and the same exponents as for t > 0. 

In the Ornstein-Zernike approximation (random-phase approximation), the 
modes with different momenta are considered to be independent of one another. 
The free energy is a sum (integral) of contributions from non-interacting normal 
modes. 



4.1.6 Validity of The Mean-field Theory 

The mean-field theories neglect the fluctuations of the neighbouring spins, each 
spin interacts only with the mean value of its neighbours. So, the mean-filed 
theories can only be valid when fluctuations are irrelevant. We will now estimate 
the relevance of fluctuations close to the critical point. We shall calculate the free 
energy of the fluctuations and compare it with the free energy calculated in the 
MFA. 

To the order of magnitude, only those fluctuations can be thermally excited 
whose free energy is of the order ksT (or less). The free energy of a typical 
fluctuation will thus be proportional to k B T. Since the size of the fluctuation is of 
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the order of the correlation length £ (this is how the correlation length is defined), 
the free energy of a typical fluctuation per unit volume is proportional to 

G fluct ex C d oc \t\ vd . (4.58) 

On the other hand, the singular part of the mean-field free energy is governed by 
the specific-heat exponent a: 

G MF oc \t\ 2 ~ a (4.59) 

(we obtain this expression after two integrations of the specific heat over temper- 
ature). For the mean-field theory to be "exact", Gf luct <C G M f in the limit as 

t -> 0, 

|t|-o G M f l*l-»o 

2 — a 

^ud-2 + a>0^d> . (4.60) 

v 

For the Ising model (scalar order parameter) in the mean-field approximation a = 
and v — 1/2 and we get the condition 

d > 4. (4.61) 

The mean-field exponents are correct and the fluctuations are irrelevant only for 
d > 4. d = 4 is called the upper critical dimension of the Ising model. 
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4.2 Monte Carlo Simulations 

With the modern computers, different computer simulations became extremely 
powerful techniques. One sets up in the computer the Hamiltonian describing a 
physical system and tries to see how it behaves. However, an efficient and reliable 
computer simulation is full of tricks and a lot of experience is needed to produce 
reliable results. 

Among different computer simulations, the Monte Carlo (MC) simulations are 
probably the most widely used in statistical mechanics. MC simulations are used 
to study classical as well as quantum systems, continuous spin systems, liquids 
and other disordered systems, gauge theories, etc. Here we shall explain the prin- 
ciples of the MC methods on the prototype example, the Ising model on a square 
lattice. 



4.2.1 Importance Sampling 

In statistical mechanics we calculate the expectation values of thermodynamic 
observables. For a spin system, the expectation value is defined as 



For an Ising system, the sum is over 2 N configurations. This number is extremely 
large for all but very small systems. Therefore one must find a method which 
estimates thermodynamic properties by sampling a small subset of representa- 
tive configurations. One possible strategy would be to scan randomly the whole 
configuration space (scan randomly over all spin configurations). The expecta- 
tion value of A would then be obtained as a sum over the scanned configurations, 
weighted by the appropriate Boltzmann factors. However, this random sampling 
has a very serious drawback. Because of the Boltzmann factor which brings a 
negligible weight to most of the configurations, very few configurations will con- 
tribute to the expectation value of A and a very unreliable estimate will result. 
This problem occurs because only an extremely restricted part of the configura- 
tion space is occupied with a considerable probability in the thermodynamic limit. 
Therefore it makes sense to restrict the sampling only to these states, this is the 
importance sampling. 

Importance sampling is realized if we generate a Markov chain of configura- 
tions, i.e., a configuration is generated from the preceding configuration. The dis- 
advantage of this Markov chain is that the successive configurations are strongly 
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correlated and therefore very long sequence of configurations is needed to obtain 
reliable averages. The expectation value of A is now equal to: 



(A) = lim 



n max n min n=n min 



(4.63) 



where A n is the average of A over n successive configurations. 



4.2.2 Metropolis Algorithm 

In generating successive configurations, the transition probability must be normal- 
ized, ergodic, and it must obey detailed balance. The Metropolis algorithm fulfills 
the above conditions. A new configuration {Sf} is obtained from the previous 
configuration {Si} by flipping one or more spins. In the Metropolis algorithm, the 
probability that the system changes from {Si} to {Sf} is 



r e -p(E f -Ei) if E > E . 

P({S t } {S,}) = 'f (4.64) 



where E { and Ef are the energies of the previous and final states, respectively. 
With this choice of transition probabilities, the system tends asymptotically to a 
steady state in which the probability of a given configuration is e^i 3E d s }\ 

The practical implementation of the Metropolis algorithm is very simple. First, 
a spin is selected either at random or sequentially. The spin is flipped and the en- 
ergies Ei and Ef axe calculated. Finally, a pseudo-random number is used to 
accept or reject the flip with the probability (4.64). The meaning of this condition 
is straightforward. If the systems gains in energy by flipping the spin, the flip 
is accepted. It is also accepted if it costs energy to flip the spin, but the energy 
difference is small so that the Boltzmann factor is larger than the random number 
(analogy to thermal excitations!). The flow chart of a MC program is shown in 
Fig. 4.6. 



4.2.3 Practical Details and Data Analysis 

• Initial Configuration. In principle, the initial configuration should not in- 
fluence the final configuration. However, the equilibrium may be reached 
faster if a proper initial configuration is chosen. An example, showing the 
approach to thermal equilibrium in a d = 3 Ising ferromagnet, is shown in 
Fig. 4.7. 
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Figure 4.6: Flow diagram of a Monte-Carlo simulation with Metropolis algorithm. 
During the first k\ steps the system still approaches to its equilibrium configura- 
tion, therefore these configurations are not included in calculating the average 
properties. 
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Figure 4.7: Relaxation of magnetization towards its equilibrium value in a Monte- 
Carlo simulation of an Ising system on a 10 x 10 x 10 lattice with periodic boundary 
conditions. Solid lines are "instantaneous" magnetizations and dashed lines their 
equilibrium estimates. (From: K. Binder and Z. Rauch, Z. Phys. 219, 201 (1969).) 
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Figure 4.8: Scaled magnetization vs. scaled temperature of an Ising model on a 
square lattice. (From: D.P. Landau, Phys. Rev. B 13, 2997 (1976).) 

• Boundary Conditions. Different boundary conditions are used. In periodic 
boundary conditions the system is periodically repeated in all directions to 
mimic an infinite system. Alternatively, one may describe the effect of the 
surrounding spins on an (inevitably) finite lattice by a mean field, acting on 
the spins at the border. Open boundary conditions, when the finite lattice is 
isolated and the spins at the surface have less neighbours than other spins, 
have also been used. 

• Finite-size effects. Lattices, considered in MC simulations, must have a 
finite size (L) whereas we are usually interested in the macroscopic prop- 
erties. The question is, how to take the limit of L to oo. Usually, one 
calculates an average for several different lattice sizes and then extrapolates 
the average to infinite L. This extrapolation is simple away from the critical 
point where the correlation length is small compared to L. The problem 
becomes non-trivial in the critical region because £ which would normally 
diverge, cannot exceed L. As a consequence, singularities associated with 
the critical point are rounded. If possible, scaling forms should be used in 
extrapolating L — >• oo (see Fig. 4.8). 

• Random-number generator. Generation of uncorrelated random numbers is 
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not a trivial problem. Computer languages have built-in "pseudo-random- 
number" generators. Not always these pseudo-random-number generators 
are random. In particular, if random sampling is used, then it might happen 
that some points on the lattice are visited more often than the other (which 

- in the worst case - may never be visited). Pseudo-random-number gen- 
erators use algorithms to calculate a sequence of (random) numbers. Care 
must be taken to initialize random number generators with different "seeds" 
in order not to obtain the same sequence of random numbers each time. 

• Statistical errors. To obtain reliable results, an average must be taken over n 
much larger than the number over which the MC states are correlated. This 
becomes more difficult near the critical point because of critical slowing 
down. 

• Free energy. Monte-Carlo simulation samples only a very small part of the 
phase space whereas to obtain the partition function or the free energy, we 
must sum over the whole configuration space. In MC simulations, there- 
fore we calculate the free energy by integrating one of its derivatives, e.g., 

— / MdH. This means that we need to calculate M(H) in a range of H . 

• Conservation laws. Above, the spins were allowed to flip uncorrelated from 
one another. This works fine for magnetic systems, where the order param- 
eter is not conserved. For binary alloys, e.g., however, the number of atoms 
is conserved and so is the order parameter. So, opposite spins have to flip 
simultaneously (one "up"^"down" and another - usually a neighbouring - 
one "down"^"up"). The conservation laws have important consequences 
for the kinetics of fluctuations (slowing down). 
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